Designing of a new multi-epitope vaccine against Leishmania major using Leish-F1 epitopes: An In-silico study

Cutaneous leishmaniasis (CL) is the most common form of the disease which can cause malignant lesions on the skin. Vaccination for the prevention and treatment of leishmaniasis can be the most effective way to combat this disease. In this study, we designed a novel multi-epitope vaccine against Leishmania major (L. major) using immunoinformatics tools to assess its efficacy in silico. Sequences of Leish-F1 protein (TSA, Leif, and LMSTI1) of L. major were taken from GenBank. The helper T (Th) and cytotoxic T (Tc) epitopes of the protein were predicted. The final multi-epitope consisted of 18 CTL epitopes joined by AAY linker. There were also nine HTL epitopes in the structure of the vaccine construct, joined by GPGPG linker. The profilin adjuvant (the toll-like receptor 11 agonist) was also added into the construct by AAY Linker. There were 613 residues in the structure of the vaccine construct. The multi-epitope vaccine candidate was stable and non-allergic. The data obtained from the binding of final multi-epitope vaccine-TLR11 residues (band lengths and weighted scores) unveiled the ligand and the receptor high score of binding affinity. Moreover, in silico assessment of the vaccine construct cloning achieved its suitable expression in E. coli host. Based on these results, the current multi-epitope vaccine prevents L. major infection in silico, while further confirmatory assessments are required.


Introduction
The worldwide outbreak of leishmaniasis, which has endangered 98 countries and a population of about 350 million people over the years, has led to many efforts to solve this health crisis [1].Unfortunately, despite many studies about leishmaniasis, there is still no definitive treatment for this zoonotic disease in human [2,3].Drugs such as pentavalent antimony (Sb 5+ ), pentamidine, amphotericin B and miltefosine are used to combat against types of the disease, but due to high cost and side effects, they cannot be prescribed with confidence [4][5][6].In addition, many vaccines were successful in vitro, but their effectiveness has not been confirmed in later stages of research and few vaccines such as LEISH-F1 and ChAd63KH DNA vaccine have entered into clinical trials [7,8].
Understanding the host defense mechanism against the parasite is essential in the development of an effective vaccine.Studies on leishmaniasis have demonstrated that provocation of cellular immunity is essential for defense against the parasite [9].Once the parasite initiates the host infection, macrophages M1 or M2 play a key role in the parasite's defense or survival.Following the T helper 1 (Th1) response elicit, the cytokines IL-12, IL-2, IL-17, IFNγ and TNF-α are increased.The presence of these cytokines, especially IFNγ, leads to the differentiation of M1 macrophages, the production of nitric oxide (NO) by inducible nitric oxide synthase (iNOS) causing the elimination of the parasite.On the other hand, by directing the immune response to T helper 2 (Th2), cytokines such as IL-4, IL-13, and IL-10 are produced, which inhibit IFN-γ production and differentiate M2 macrophages.This mechanism causes the infection persistence in the body due to the inefficiency of the immune system in parasitic defense [10,11].Also, in some studies, an increase in humoral immune response has been demonstrated by IgG2a antibody production to defense against leishmaniasis [12].Therefore, vaccination aimed at stimulating Th1, Tc, and humoral immune responses can be highly effective in defense against the parasite.
Several Leishmania antigens have been used to synthesize a variety of vaccines against the disease, some of which have inferred sufficient immunogenicity in vitro.One of these antigens is Leish-F1 containing Leishmania major (L.major) Thiol-specific Antioxidant (TSA), L. major stress inducible 1 (LMSTI1), and L. braziliensis elongation initiation factor (Leif) [13].The immunogenicity evaluation has been performed using a vaccine containing Leish-F1 antigen with MPL-SE adjuvant and all of them have been successful in clinical trials [14].
These immunogens can be used to make a variety of vaccines, including protein, DNA and mRNA vaccines.As mentioned, the Leish-F1 protein has demonstrated effectiveness with potential of the immune system provocation, making it a prominent candidate for synthesizing other types of vaccines.On the other hand, due to the large size of the genome of this polyprotein, adoption of potential epitopes is reliable.Considering as a candidate for the synthesis of the DNA vaccine against Leishmania spp, it can be cloned into the AAV (adeno-associated virus) vectors [15].
Nowadays, with the advancement of immunoinformatics approach, vaccines can be designed in silico prior to in vitro and in vivo development for the efficacy assessment.According to the prediction methods of immunoinformatics, the most suitable epitopes of each immunogens are selected for vaccine design.Thus, it can decrease the cost and duration of vaccine development along with increase of safety, specify, and effectiveness of vaccines [16,17].Also, using a multi-epitope instead of a whole protein vaccine reduce allergic reactions enhancing a protein vaccine efficacy [18,19].Moreover, due to the lack of living infectious agents in the design these multi-epitope vaccines [20,21], the risk of infection among vulnerable individuals is substantially mitigated in this preventive approach [22].
In this study, for the first time, we designed a multi-epitope vaccine using Leish-F1 immunogenic antigen of L. major to combat against the CL using immunoinformatics tools.To increase the immunogenicity of the vaccine, we also used the profilin adjuvant, the toll-like receptor 11 (TLR11)-agonist, which can interact with the macrophages and monocytes enhancing their recruitment and activation.For this purpose, the most suitable HTL and CTL epitopes of Leish-F1 were predicted and linked by AYY and GPGPG linkers.Also, Profilin with the AAY linker was added to the selected epitopes and predicted in terms of three dimensional (3D) structure validity and refinement, and B-cell epitopes prediction.Furthermore, molecular docking of final vaccine and TLR11, molecular dynamics (MD) simulation and codon optimization were performed.

T-cell epitope prediction
To predict the CTL and HTL epitopes of TSA, LMSTI1, and Leif proteins, the IEDB server was used.Prediction of 9-mer and 15-mer epitopes to bind to MHCI and MHCII were done by using IEDB (http://tools.iedb.org/mhcii/).The server predicts the CTL and HTL epitopes on the basis of SMM (stabilized matrix method), Artificial Neural Networks (ANN) and, Allele specific affinity cutoff (IC50) [23].In this study, we selected the mouse alleles to bind the MHCI and MHCII.

IFN-γ inducing epitope prediction
Using IFN epitope server (http://crdd.osdd.net/raghava/ifnepitope/),prediction of gammainterferon inducing epitopes of HTL epitopes was performed.The server predicts the epitopes based on the methods such as SVM based, motif based, and hybrid approach.On this server, IFN-γ inducing and non-inducing MHC-II binders were determined.

Construction of final multi-epitope vaccine
Using the in silico information, the final multi-epitope vaccine was designed.First, AYY and GPGPS flexible linkers were used to join HTL and CTL epitopes, respectively [24].Then, to enhance the immunogenicity of final multi-epitope vaccine, profilin was added as an adjuvant to the N-terminal of the construct using AYY linker.

Conformational and linear B cell epitope prediction
According to the 3-D structure of the construct, the ElliPro suite (http://tools.iedb.org/ellipro/ ) was used to predict the conformational and linear B cell epitopes of the final vaccine construct [25].The server assigns a score to each predicted epitope, named PI (Protrusion index).Firstly, the three-dimensional structure of the protein is estimated by a number of spherical ellipses, which the higher the percentage, the fewer residues outside the ellipsoid.The PI value is defined based on the center of mass remaining outside the largest possible ellipsoid and residues with higher scores have more access to the solvent.Discontinuous epitopes are characterized by PI values and clustered by R distance (in Å between remaining mass centers).The higher the R value, the larger the epitopes.

Prediction and validation of protein structures of TLR11 and final multi-epitope vaccine
AlphaFold2 program was used to predict of tertiary structure of the TLR11 and final vaccine construct, based on a deep neural network algorithm [26].PROCHECK (http://servicesn.mbi.ucla.edu/PROCHECK/),ProSA (http://prosa.services.came.sbg.ac.at/prosa.php),and ERRAT (http://services.mbi.ucla.edu/ERRAT)servers were used to validate the 3D structure of the TLR11 and final vaccine construct [27].The PROCHECK gave a Ramachandran plot and assessed the stereochemical quality of the refined model [28].Furthermore, the ProSA provided an overall quality score based on Z-score value.It should be noted that scores outside the range of native proteins are not acceptable [29].

Evaluation of physico-chemical properties of the final vaccine construct
ProtParam (http://web.expasy.org/protparam/),an online web server was used to calculate the physicochemical characteristics of final multi-epitope vaccine containing aliphatic index (AI), instability index, theoretical PI (isoelectric point) molecular weight, in vitro and in vivo halflife, and grand average of hydropathicity (GRAVY) [30].Then, the solubility of the final vaccine construct over expression in E. coli was calculated using SOLpro (http://scratch.proteomics.ics.uci.edu/)[31].

Allergenicity and antigenicity assessment of the final vaccine construct
AllerTOP v2.0 and AllergenFP1.0were applied to determine the antigenicity of the final vaccine construct.Using five E-descriptors, including amino acid hydrophobicity, size, amino acid tendency to helix, amino acid abundance, and tendency of β-strand formation, the amino acids in the protein sequence are described in the data set of AllergenFP1.0server.Also, based on auto-cross covariance (ACC) transformation, the strings were transformed into uniform vectors.The server categories the proteins in two allergen and non-allergen, according to Tanimoto coefficient [32].AllerTOP v2.0 also used the same five aforementioned E-descriptors.This server classifies the protein to allergen and non-allergen based on k-nearest neighbor algorithm (kNN, k = 1) [33].Thereafter, the determination of antigenicity was done using ANTIGENpro (http://scratch.proteomics.ics.uci.edu/) and VaxiJen v.2 (http://www.ddgpharmfac.net)servers.ANTIGENpro is an alignment-free and sequence-based server.The accuracy of this pathogen independent server is 82% [34].Protective antigens were also predicted by VaxiJen, a first alignment-independent prediction [35].

Docking analysis of vaccine construct and TLR11
ClusPro server available at (http://cluspro.bu.edu/login.php) was applied to determine the affinity of binding between ligand (final multi-epitope vaccine) and receptor (TLR11).Elimination of unstructured protein regions, use of attraction or repulsion, calculation of pairwise distance restraints, and homo-multimers construction, are some of the parameters that this server considers to provide docking results of two proteins.In addition, the selection of highly populated cluster centers of low-energy structures is another feature of this server [36].Finally, visualization of docked complex was performed by PyMOL program.

Molecular dynamics simulation details
The structural properties of the TLR11, multi-epitope vaccine, the final vaccine construct-TLR11 and the TLR11-adjuvant (as positive control) complexes were evaluated in the in silico conditions through MD simulation.In this step, Using GROMOS 54a7 and 43a1 package, docked complex of final vaccine and TLR-11 was simulated during 50000 picoseconds (Ps) at 310˚K and 1 bar pressure.A Cubic box with periodic boundary conditions with water molecules of SPC/E was also used.Moreover, to minimize the system energy and relaxation of solvent molecule, the steepest-descent algorithm was applied.Also, fixation of the atom was performed by using A LINear Constraint Solver (LINCS) algorithm and SETTLE, an analytical algorithm, was used for solvent molecules.To calculate of total electrostatic energy and other non-bonded interactions, Particle Mesh Ewald (PME) summation method and L-J model (with 10 A˚cutoff distances), were used, respectively.Moreover, the temperature and pressure of each system, at the time of coupling equal to 0.1 Ps were preserved using Berendsen weak coupling algorithm.Finally, the root mean square deviation (RMSD), root mean square fluctuation (RMSF), hydrogen bond (H-bond), and radius gyration (Rg) diagrams were plotted for both systems.

In silico cloning of final multi-epitope vaccine
For this purpose, SnapGene tool was applied and cDNA was provided as an output of the server.The cDNA sequence was used to determine the CAI value and GC content.The CAI value higher than 0.8 and GC content between 30-70% are acceptable.Then, XhoI and NdeI restriction enzymes were added to the cDNA sequence.Finally, codon optimization into the pET28a (+) vector was done.

T-cell epitope prediction
Using IEDB server, CTL epitopes of leishF1consisting three proteins was predicted.Six epitopes were predicted from each protein of leishF1.Therefore, a total of 18 epitopes with the lowest percentile rank was selected.Also, a total of nine HTL epitopes of all proteins was selected and for each protein three epitopes were selected (Tables 1 and 2).Moreover, applying IFNepitope server, 12 HTL epitopes were determined to induce of IFN-γ.(Table 3).

Construction of final multi-epitope vaccine
A total of 18 CTL epitopes was linked together by AAY linker.Also, GPGPG linker was utilized to join of 9 HTL epitopes with them.Then, profilin with the length of 567 amino acid was added at the N terminal of the construct.The construct was shown in Fig 1.

B cell epitope prediction of final multi-epitope vaccine
Prediction of linear and discontinuous epitopes of the final vaccine construct was performed using ElliPro suite.According to the results obtained from the server, 9 linear and 11 discontinuous B cell epitopes were predicted in final vaccine with the potential of stimulating humoral immunity (Fig 2A and 2B and Tables 4 and 5).

Prediction and validation of three-dimensional structure of TLR11 and final multi-epitope vaccine
AlphaFold2 results showed five predicted models for the 3D structure of the TLR11 and vaccine.These predicted structures were validated using methods such as the Ramachandran plot, MolProbity, ProSA, and ERRAT.For this purpose, verification of quality and potential errors of the crude 3D models was done.ERRAT method showed values of 81.25% and 80.40% as an overall quality factor of the best models (3A and B), for TLR11 and vaccine respectively (Table 6).The ProSA calculated a Z-score of -4.46, as a plausible value of native proteins as   this step such as half-life in mammalian reticulocytes were 30 hours in mammalian reticulocytes, in vitro, >20 hours in yeast and >10 hours in E. coli and in vivo.Also, instability index and AI were computed 29.17 and 75.22,respectively.The stability and thermostability of the construct in nature was confirmed by the values of attained from instability index and AI.The GRAVY score of -0.191 unveiled the construct hydrophilicity in nature.Then, SOLpro server was utilized to investigate the final vaccine construct solubility and antigenicity.The solubility of the vaccine candidate in E. coli and its antigenicity were confirmed giving values of 0.91 and 0.72 respectively.

Evaluation of allergenicity and antigenicity of final multi-epitope vaccine
Allergenicity assessment was done using AllergenFP1.0and AllerTOP servers.The results revealed the final vaccine construct non-allergenicity.Moreover, its antigenicity was carried

Docking study between final multi-epitope vaccine and TLR-11
First of all, the MD simulations of the final multi-epitope vaccine and TLR-11 predicted models were performed for 50000 ps, then, ClusPro v. 2.0 webserver was used for the generation of the complex systems.Based on biophysical properties of ligand (multi-epitope vaccine and adjuvant) and receptor (profilin) a total of 30 models were predicted.Among them, the model number 0.00 with the lowest energy score was selected.The weighted score obtained were determined -801.1 (vaccine-TLR11) and -711.5 (adjuvant-TLR11) kcal with the highest binding affinity in all of the predicted models.These complexes were used as initial structures for MD simulation.After 50 ns MD simulation, for a better understanding of the protein-protein interface region of the complex systems, pdbsum web server was used.These analyses show that there are hydrogen bonds, salt bridges, and non-bonded contacts at the interface region of the vaccine-TLR11 and adjuvant-TLR11 complexes (Fig 4A and 4B).Based on this analysis, 9 salt bridges, 21 hydrogen bonds, and 173 non-bonded contacts were found in the interface region of the vaccine-TLR11 complex.Although, 5 salt bridges, 11 hydrogen bonds and 85 non-bonded contacts were found between adjuvant and TLR11 (Fig 4A and 4B).

Molecular dynamics simulation trajectory analysis
The RMSD values for the multi-epitope vaccine and TLR11 in the free state were assessed as a time-dependent parameter for calculating the displacement of α-carbon atoms in the two molecular types during the MD simulations (Fig 5A).In accordance with the RMSD diagrams, the TLR11 has a small fluctuation compared to that of the multi-epitope vaccine.The mean RMSD values for the TLR11 and multi-epitope vaccine were respectively (0.22 ± 0.01) nm and (1.3 ± 0.3) nm.Also, the RMSF values of the TLR11 and the multi-epitope vaccine were estimated separately for analyzing the local structural fluctuations (Fig 5B).High values of the RMSF diagrams refer to loop regions in the protein structures in two proteins.The compactness of the protein structures (TLR11 and vaccine) was estimated using the Rg method.The Rg graph of the vaccine (Fig 5C ) presented that over the simulation period from 0 ns to 10 ns with sharp slope decreases from 3.4 nm to 3.0 nm, after that, the graph with a slight slope decreased.This indicated that the multi-epitope vaccine changes its configuration more rapidly in this period.It also can be seen in Fig 5A that the sharp increase of RMSD is in accordance with the curve of Rg when the compactness of the multi-epitope vaccine is low.With the increase of protein compactness, the RMSD graph reaches a plateau.Our results unraveled that the multi-epitope vaccine structure compactness during a simulation time increases which in turn reflects its stability.Rg and RMSD fluctuations of the TLR11 are very tiny which presented the high structural stability.Also, Gibb's free energy was analyzed to study the conformational stability of the modeled vaccine using the first two eigenvectors (PC1 vs PC2) from the simulation trajectory.According to Fig 5D, the energy conformations of the modeled vaccine were limited to the narrow range of energy basins -40 to 10 kJ/mol on PC1 and -20 to 10 kJ/mol on PC2.The RMSD value between the model protein structure and spatial configuration with the least energy of the vaccine structure was calculated to be 1.06 nm (Fig 5E), which may be attributed to the comparatively lower compactness of the model protein structure in comparison to the spatial configuration exhibiting the minimum energy in the vaccine structure.Alternatively, the presence of loop regions in both structures could also lead to this Our results outlined that compared to the adjvant-TLR11 as positive control, the vaccine complex had stable dynamics with the TLR11 receptor, confirming the strong binding of the vaccine to the TLR11.

Contribution of energy components to the vaccine binding along MD simulations
The interaction energies between each of the vaccine and adjuvant with TLR11 was estimated as separate electrostatic (ΔE ele ) and van der Waals (ΔE vdW ) components (Table 7) using MMPBSA approach.According to Table 7, hydrophobic (ΔE non-polar ) energy term has a favorable role in the binding each of the vaccine and adjuvant to the TLR11.Based on our results, non-polar component plays the major role in the vaccine and adjuvant binding as main driving force (Table 7).Notably, the ΔG bind of the vaccine-TLR11 complex was considerably better than that of adjuvant-TLR11, highlighting stronger and more stable binding and higher affinity of vaccine-TLR11 complex.According to the Fig 5K, the number of H-bonds between the vaccine and TLR11 (Fig 4A and 4B) in the interface regions was higher than those of adjuvant and TLR11 during the simulation time, which also confirmed the complex stability.

In silico codon optimization of final multi-epitope vaccine
Based on the results of the Java Codon Adaptation tool (JCat), reverse translation and codon optimization in E. coli (stain K12) was carried out.After the design of cDNA, codon adaptation index and GC percent were calculated 1.0 and 53.35, respectively.The value of GC content was in acceptable rang of the results, being 30-70%.Also the CAI value showed the vaccine highest expression level in the pET28a (+) vector (Fig 6).

Discussion
There is no definitive treatment for leishmaniasis, in particular the CL form of the disease, which affects 1.5-2 million people annually, and can cause several lesions on the patient's skin if it progresses, may even lead to the formation of bacterial biofilms and antibiotic resistance in the patient [37][38][39].Researchers attempt to make a vaccine for prevention and elimination of the disease [40].Identifying the host defense pathways against the pathogen and then strengthening these pathways can be the most important strategy in making an effective vaccine.To achieve this goal, immunoinformatics tools can be used to select the most appropriate epitope of the desired immunogens in less time with high accuracy [17,41].
The aim of this study was to design a novel multi-epitope vaccine to be effective in combating leishmaniasis.In the first step, HTL and CTL epitopes of Leish-F1 polyprotein from L. major were predicted and ligated using AYY and GPGPG linkers.Also, due to the fact that in the vaccination against leishmaniasis, the immune system needs to be sufficiently activated, we used the profilin as an adjuvant.The Profilin can bind to TLR11 at the macrophage and dendritic cells surface, leading to innate immune stimulation, increased IFN-γ, and consequently increased Th1.Therefore, profilin can help increase the immunogenicity of this vaccine against L. major, which was our main goal [48].There were also 12 epitopes in the final vaccine construct to induce INF-γ, an indispensable cytokine eliciting the immune system response against Leishmania spp, these epitopes can be promising for the defense against the parasite.
After the design of the final vaccine construct, conformational and linear B cell epitopes were predicted.The number of twenty linear and conformational epitopes in the structure of the final construct showed that this construct can properly provoke the humoral immunity of the host.
Moreover, tertiary structure prediction of final vaccine construct along with its validation and refinement was implemented.Based on the results of Ramachandran diagram and ProSA, the final vaccine construct was stable in nature.Also the quality and potential errors were predicted using the refinement of final vaccine.
Thereafter, on the basis of physicochemical properties of final vaccine, the PI and GRAVY values included 5.51 and -0.191, respectively, indicating the construct acidity and hydrophilicity in nature.Also, thermostability of the final vaccine was confirmed with the values of 75.22 and 29.17 for AI and instability index, respectively.Then, the solubility of the final multi-epitope vaccine in E. coli was confirmed.
According to the Rg, RMSD, RMSF and H-bond values obtained from MD trajectory both of the vaccine-TLR11 and adjvant-TLR11 complexes had well structural stability.Additionally, MMPBSA results demonstrated that the hydrophobic residues play a main role in the binding process between each of the vaccine and adjuvant with the TLR11, which in turn maintains the high stability of the complex structures.Therefore, the findings revealed acceptable affinity of the designed multi-epitope vaccine to the TLR11, exhibiting high complex structural stability compared to the adjvant-TLR11 as positive control.The size of the vaccine construct is higher that adjuvant, therefore, more residues interact to the TLR11, and the difference in the interface residues is the result of vaccine and adjuvant differences.
In many vaccine design studies, docking and MD simulations are performed to ensure the efficacy of the construct.Because by confirming the affinity and stability of ligand-receptor binding in conditions similar to the body environment, the efficacy of the vaccine in the in silico conditions is confirmed.Finally, the results of codon optimization (CAI value and GC content) showed that the multi-epitope vaccine could have the highest expression in the host.

Conclusion
In this study, using immunoinformatics tools, a novel multi-epitope vaccine was designed consisting the Leish-F1 of L. major.The vaccine construct was designed using the most suitable epitopes of cellular immune stimulation of the host and adding an adjuvant (TLR11 agonist) to increase its immunogenicity.The B cell results showed that there are epitopes of humoral immune stimulation in the final construct since humoral immunity can be also effective in defense against leishmaniasis.Moreover, molecular docking and the MD simulation analysis of final multi-epitope vaccine and TLR11 compared to the adjvant-TLR11 as positive control highlighted acceptable binding free energy and confirmed the stability of the vaccine-TLR11 complex.Based on these results, the multi-epitope vaccine was non-toxic and non-allergenic, immunogenic, antigenic, soluble and stable and thereby could stimulate the immune system against L. major.Despite the vaccine efficiency in the in silico, it is necessary to confirm its effectiveness in the further studies in vitro and in vivo.
shown in Fig 3C.In addition, the data obtained from the Ramachandran plot of the vaccine exhibited 70.7% residues located in most favored regions, 26.1% residues in allowed regions, and 3.2% in disallowed regions as shown in Fig 3D.The Z-score included -8 for the TLR11 (Fig 3E).Also, in the refined model of the TLR11, 81.1% of the residues were in the favored area, 17.9% in the allowed area and 1.1% in the disallowed area of Ramachandran plot (Fig 3F).

Fig 3 .
Fig 3. Prediction and validation of the 3D structure of the TLR11 and multi-epitope vaccine; A and B Final predicted models of the vaccine and TLR11 respectively; C and D ProSA and Ramachandran plots for the vaccine predicted model respectively; E and F ProSA and Ramachandran plots for the TLR11 predicted model respectively.https://doi.org/10.1371/journal.pone.0295495.g003